%% This file is called to compute comparative static results; 
%%%% It needs "par" as parameters
%%%% It outputs results stored in variables named similar to Inv_model_5
%%%% below

disp('-------------> Model 5, idiosyncratic shocks + imperfect reallocation + imperfect credit <---------')

%% Initialization
Inv_model_5                 = zeros(length(ALPHA_grid), length(THETA_grid));
Prod_model_5                = Inv_model_5;
Prod_std_model_5            = Inv_model_5;
Output_model_5              = Inv_model_5;
Hours_model_5               = Inv_model_5;
R_ratio_model_5             = Inv_model_5;
P_share_model_5             = Inv_model_5;
Aqc_model_5                 = Inv_model_5;
Aqc_P_model_5               = Inv_model_5;
Sppe_model_5                = Inv_model_5;
Sppe_P_model_5              = Inv_model_5;
Util_model_5                = Inv_model_5;
Welfare_model_5             = Inv_model_5;
Consumption_model_5         = Inv_model_5;
Sec_mkt_price_model_5       = Inv_model_5;
Z_Y_model_5                 = Inv_model_5;
Z_K_model_5                 = Inv_model_5;
debt_model_5                = Inv_model_5;
Threshold_model_5           = Inv_model_5;
perfect_reallocation        = 0;
cali_step                   = 0;

%% Comparative statics
for ii = 1 : length(THETA_grid)
    par = par_ini;
    par.THETA       = THETA_grid(ii);
    fun             = define_function_fun(par);
    
    fprintf('I am doing the case of THETA = %.3f\n\n', THETA_grid(ii))
    
    Inv_temp            = zeros(length(ALPHA_grid), 1);
    Prod_temp           = Inv_temp;
    Prod_std_temp       = Inv_temp;
    Output_temp         = Inv_temp;
    Hours_temp          = Inv_temp;
    R_ratio_temp        = Inv_temp;
    P_share_temp        = Inv_temp;
    Aqc_temp            = Inv_temp;
    Aqc_P_temp          = Inv_temp;
    Sppe_temp           = Inv_temp;
    Sppe_P_temp         = Inv_temp;
    Util_temp           = Inv_temp;
    Welfare_temp        = Inv_temp;
    Consumption_temp    = Inv_temp;
    Sec_mkt_price_temp  = Inv_temp;
    Z_Y_temp            = Inv_temp;
    Z_K_temp            = Inv_temp;
    debt_temp           = Inv_temp;
    Threshold_temp      = Inv_temp;
    
    parfor jj = 1: length(ALPHA_grid)
        par_temp = par;
        par_temp.ALPHA  = ALPHA_grid(jj);%
        
        x       = fsolve(@(x)eqm_ss_iid_cases_2_to_5(x, par_temp, fun, perfect_reallocation, cali_step), [0.05, par_temp.ALPHA_0], options);
        B       = x(1);
        ALPHA_0 = x(2);
        
        [fval, eqm]       = eqm_ss_iid_cases_2_to_5(x, par_temp, fun, perfect_reallocation, cali_step);
        eqm.IM            = average_log_normal(eqm.surplus_net, par.c_mu, par.c_sig) / par.XI_IM;

        Consumption_temp(jj)   = eqm.c_star + eqm.w * eqm.IM * par_temp.ZETA;
        Inv_temp(jj)      = eqm.K * par_temp.DELTA; %no normalization
        Hours_temp(jj)    = eqm.H;

        Output_temp(jj)   = (eqm.Y + eqm.w * eqm.IM * par.ZETA) / Y_Hansen;
        
        Prod_temp(jj)     = eqm.K_bar_K / K_bar_K_Hansen;
        Prod_std_temp(jj) = eqm.prod_std / prod_std_Hansen;
        R_ratio_temp(jj)  = (eqm.E_a_K + eqm.E_p_K) / (eqm.E_a_K + eqm.E_p_K + par_temp.DELTA);
        P_share_temp(jj)  = eqm.E_p_K / (eqm.E_a_K + eqm.E_p_K);
        Aqc_temp(jj)      = eqm.E_a_K  / (eqm.E_a_K + eqm.E_p_K + par_temp.DELTA);
        Aqc_P_temp(jj)    = eqm.prob_A;
        Sppe_temp(jj)     = eqm.E_p_K  / (eqm.E_a_K + eqm.E_p_K + par_temp.DELTA);
        Sppe_P_temp(jj)   = eqm.prob_P;
        Util_temp(jj)     = log(eqm.c_star) - par_temp.XI * eqm.H;
        %Welfare_temp(jj)  = eqm.c_star / c_star_Hansen * exp(par_temp.XI * (H_Hansen - eqm.H)) - 1;
        Welfare_temp(jj)  = eqm.c_star / eqm_calibrated.c_star * exp(par_temp.XI * (eqm_calibrated.H - eqm.H) + par_temp.XI_IM * (eqm_calibrated.IM - eqm.IM)) - 1; % use calibrated c_star and H

        Sec_mkt_price_temp(jj) = eqm.average_2nd_price;
        Z_Y_temp(jj)      = ALPHA_0 * eqm.Z / eqm.Y;
        Z_K_temp(jj)      = ALPHA_0 * eqm.Z / eqm.K;
        debt_temp(jj)     = ALPHA_0 * par_temp.ALPHA * eqm.debt / eqm.Y;
        Threshold_temp(jj)= eqm.threshold_iota;        
 
    end
    Inv_model_5(:,ii)      = Inv_temp;
    Prod_model_5(:,ii)     = Prod_temp;
    Prod_std_model_5(:,ii) = Prod_std_temp;
    Output_model_5(:,ii)   = Output_temp;
    Hours_model_5(:,ii)    = Hours_temp;
    R_ratio_model_5(:,ii)  = R_ratio_temp;
    P_share_model_5(:,ii)  = P_share_temp;
    Aqc_model_5(:,ii)      = Aqc_temp;
    Aqc_P_model_5(:,ii)    = Aqc_P_temp;
    Sppe_model_5(:,ii)     = Sppe_temp;
    Sppe_P_model_5(:,ii)   = Sppe_P_temp;
    Util_model_5(:,ii)     = Util_temp;
    Welfare_model_5(:,ii)  = Welfare_temp;
    Z_Y_model_5(:,ii)      = Z_Y_temp;
    Z_K_model_5(:,ii)      = Z_K_temp;
    debt_model_5(:,ii)     = debt_temp;
    Threshold_model_5(:,ii)= Threshold_temp;
    Consumption_model_5(:,ii)   = Consumption_temp;
    Sec_mkt_price_model_5(:,ii) = Sec_mkt_price_temp;    
end